use ${newdata}carownership_dataset_bergen, clear

/* Figure 3: yvar is pr(BEV) */
global yvar bev

/* Only keep required obs */
keep if $trmgroup != 0 
keep if $trmgroup != .
keep if 	(year >= $firstyearpre  & year <= $lastyearpre) ///
		  | (year >= $firstyearpost & year <= $lastyearpost)
		  
/* setup variables */
capt drop bergen 		  
gen commuters = ($trmgroup == 1 | $trmgroup == 3) 
gen bergen = 	($trmgroup == 1 | $trmgroup == 2)
gen post = (year >= $firstyearpost & year <= $lastyearpost)
gen bergen_commuters = bergen * commuters
gen commuters_post = commuters * post
gen bergen_commuters_post = bergen_commuters * post



/* All time-varying variables based on geography (residence or work location)
   are set to their 2014 values (variation post 2014 is endogenous to treatment) */
foreach var in dist time_work PublicVSCarTime_fam_mean PublicDiffCarTime_fam_mean grk {
	replace `var' = . if year != 2014
	bysort familienr (`var'): replace `var' = `var'[_n-1] if missing(`var')
}	  

sort familienr year
egen long famid_num = group(familienr)

/* Main work type */
gen main_educ1 = substr(yrk_kode1,1,1)
destring main_educ1, replace
gen main_educ2 = substr(yrk_kode2,1,1)
destring main_educ2, replace
	
/* === First heterogeneous treatment: Household income ====================== */		  
gen heterovar_inc = 0
replace heterovar_inc = wies1 if wies1 != .
replace heterovar_inc = heterovar_inc + wies2 if wies2 != .	
 
xtile heterovar_inc_qt = heterovar_inc if year == 2014, nquantiles(5)
/*
table heterovar_inc_qt, stat(mean heterovar_inc)
table heterovar_inc_qt, stat(min heterovar_inc)
table heterovar_inc_qt, stat(max heterovar_inc)
*/

capt drop aux
bys famid_num: egen aux = mean(heterovar_inc_qt)
replace heterovar_inc_qt = aux
drop aux
	
forvalues i = 1/5 {
	gen bergen_commuters_post_`i' = bergen_commuters_post * (heterovar_inc_qt == `i')
}			  

/* Regression */
eststo clear
reghdfe $yvar ///
	bergen_commuters_post_* /// DiDiD: post commuter in bergen
	i.year##i.heterovar_inc_qt /// 
	$xvar ///
	i.commuters##i.heterovar_inc_qt##i.year /// 
	i.bergen_commuters##i.heterovar_inc_qt ///
	i.bergen##i.heterovar_inc_qt ///
	, absorb( ///
		i.grk#i.year ///
		i.famid_num ///
	) vce(cluster grk)	
eststo reg_inc
forvalues i = 1(1)5 {
	summarize heterovar_inc if e(sample) == 1  & heterovar_inc_qt ==`i'
	estadd scalar MeanHetero_`i' = r(mean), replace
	estadd scalar NHetero_`i' = r(N), replace
	summarize $yvar if e(sample) == 1  & year == $lastyearpre ///
	& bergen_commuters == 1 & heterovar_inc_qt ==`i'
	estadd scalar MeanValue1_`i' = r(mean), replace
	summarize $yvar if e(sample) == 1  & year == $lastyearpost ///
	& bergen_commuters == 1 & heterovar_inc_qt ==`i'
	estadd scalar MeanValue2_`i' = r(mean), replace
}
eststo reg_inc 
estimates save ${ster}reg_hetero_${yvar}_income, replace

/* Loading estimates manually */
estimates use ${ster}reg_hetero_${yvar}_income

/* === FIGURE of the last of the regressions === */
/* Saving estimates */
matrix B = e(b)
matrix V = e(V)

/* Storing them as variables */
capt drop xaxis beta se meanheterovar 
gen xaxis = .
gen beta = .
gen se = .
gen meanheterovar = .

forvalues i = 1/5 {
	qui replace xaxis = `i' in `i'
	qui replace beta  = B[1,`i'] in `i'
	qui replace se    = sqrt(V[`i',`i']) in `i'
	qui sum heterovar_inc if heterovar_inc_qt == `i' & year == $lastyearpre ///
	& bergen_commuters == 1
	qui replace meanheterovar = r(mean) in `i'
}

/* 95 % CI */
capt drop plusse minusse
qui gen plusse = beta + 1.96 * se
qui gen minusse = beta - 1.96 * se
replace meanheterovar = meanheterovar / 1000

/* === FIGURE ============================================================*/

/* Table to find correct labels for figure */
table heterovar_inc_qt if year == 2014, stat(min heterovar_inc) stat(max heterovar_inc)

twoway ///	
	(rcap plusse minusse xaxis ///
		, sort color(black)) ///			
	(scatter beta xaxis ///
		, sort lcolor(black) mcolor(black) lpattern(solid)) ///	
	, graphregion(color(white))	///
	xtitle("Income quintile (1000 NOK)") ytitle("Treatment effect") ///
	name(fig3a, replace) ///
	xlabel(1 "< 393" 2 "393-575" 3 "575-727" 4 "727-905" 5 ">905", angle(55)) ///
	legend(off) yline(0) ///
	ylabel(-0.02(0.02)0.08) ///
	scale(1.3)

graph export "${figures}fig5a.png", as(png) replace
graph export "${figures}fig5a.pdf", as(pdf) replace

/* === Second heterogeneous treatment: age ================================== */		  

  
xtile heterovar_age_qt = age if year == 2014, nquantiles(5)
*table heterovar_age_qt, stat(max age)
bys famid_num: egen aux = mean(heterovar_age_qt)
replace heterovar_age_qt = aux
drop aux

drop bergen_commuters_post_*	
forvalues i = 1/5 {
	gen bergen_commuters_post_`i' = bergen_commuters_post * (heterovar_age_qt == `i')
}			  



eststo clear
reghdfe $yvar ///
	bergen_commuters_post_* /// DiDiD: post commuter in bergen
	i.year##i.heterovar_age_qt /// 
	$xvar ///
	i.commuters##i.heterovar_age_qt##i.year /// 
	i.bergen_commuters##i.heterovar_age_qt /// 
	i.bergen##i.heterovar_age_qt ///
	, absorb( ///
		i.grk#i.year ///
		i.famid_num ///
	) vce(cluster grk)	
eststo reg_3
forvalues i = 1(1)5 {
	summarize age if e(sample) == 1  & heterovar_age_qt ==`i'
	estadd scalar MeanHetero_`i' = r(mean), replace
	estadd scalar NHetero_`i' = r(N), replace
	summarize $yvar if e(sample) == 1  & year == $lastyearpre ///
	& bergen_commuters == 1 & heterovar_age_qt ==`i'
	estadd scalar MeanValue1_`i' = r(mean), replace
	summarize $yvar if e(sample) == 1  & year == $lastyearpost ///
	& bergen_commuters == 1 & heterovar_age_qt ==`i'
	estadd scalar MeanValue2_`i' = r(mean), replace
}
eststo reg_3 
estimates save ${ster}reg_hetero_${yvar}_age, replace
/* === FIGURE of the last of the regressions === */
/* Saving estimates */
matrix B = e(b)
matrix V = e(V)

/* Storing them as variables */
capt drop xaxis beta se meanheterovar 
gen xaxis = .
gen beta = .
gen se = .
gen meanheterovar = .
*local coef = $fig4lastyear - $fig4firstyear /* The number of coef we want to keep */
forvalues i = 1/5 {
	qui replace xaxis = `i' in `i'
	qui replace beta  = B[1,`i'] in `i'
	qui replace se    = sqrt(V[`i',`i']) in `i'
	qui sum age if heterovar_age_qt == `i' & year == $lastyearpre ///
	& bergen_commuters == 1
	qui replace meanheterovar = r(mean) in `i'
}

/* 95 % CI */
drop plusse minusse
qui gen plusse = beta + 1.96 * se
qui gen minusse = beta - 1.96 * se

/* === FIGURE ============================================================*/
table heterovar_age_qt if year == 2014, stat(min heterovar_inc) stat(max heterovar_inc)
	
twoway ///	
	(rcap plusse minusse xaxis ///
		, sort color(black)) ///			
	(scatter beta xaxis ///
		, sort lcolor(black) mcolor(black) lpattern(solid)) ///	
	, graphregion(color(white))	///
	xtitle("Age quintile") ytitle("Treatment effect") ///
	name(fig3d, replace) ///
	xlabel(1 "18-32" 2 "32-39" 3 "39-46" 4 "46-55" 5 "55-90", angle(55)) ///
	legend(off) yline(0) ///
	ylabel(-0.02(0.02)0.08) /// ylabel(-0.04(0.04)0.12) ///
	scale(1.3)	


graph export "${figures}fig5d.png", as(png) replace
graph export "${figures}fig5d.pdf", as(pdf) replace

/* === Third heterogeneous treatment: family status ========================= */		  

gen famstatus = .
replace famstatus = 1 if couple == 0 & children == 0 & year == 2014
replace famstatus = 2 if couple == 0 & children == 1 & year == 2014
replace famstatus = 3 if couple == 1 & children == 0 & year == 2014
replace famstatus = 4 if couple == 1 & children == 1 & year == 2014 
bys famid_num: egen heterovar_famstatus = mean(famstatus)

drop bergen_commuters_post_*	
forvalues i = 1/4 {
	gen bergen_commuters_post_`i' = bergen_commuters_post * (heterovar_famstatus == `i')
}			  

eststo clear
reghdfe $yvar ///
	bergen_commuters_post_* /// DiDiD: post commuter in bergen
	i.year##i.heterovar_famstatus /// 
	$xvar ///
	i.commuters##i.heterovar_famstatus##i.year /// 
	i.bergen_commuters##i.heterovar_famstatus /// 
	i.bergen##i.heterovar_famstatus  ///
	, absorb( ///
		i.grk#i.year ///
		i.famid_num ///
	) vce(cluster grk)	
eststo reg_3
forvalues i = 1(1)4 {
	summarize heterovar_famstatus if e(sample) == 1  & heterovar_famstatus ==`i'
	estadd scalar MeanHetero_`i' = r(mean), replace
	estadd scalar NHetero_`i' = r(N), replace
	summarize $yvar if e(sample) == 1  & year == $lastyearpre ///
	& bergen_commuters == 1 & heterovar_famstatus ==`i'
	estadd scalar MeanValue1_`i' = r(mean), replace
	summarize $yvar if e(sample) == 1  & year == $lastyearpost ///
	& bergen_commuters == 1 & heterovar_famstatus ==`i'
	estadd scalar MeanValue2_`i' = r(mean), replace
}
eststo reg_3 
estimates save ${ster}reg_hetero_${yvar}_famstatus, replace
/* === FIGURE of the last of the regressions === */
/* Saving estimates */
matrix B = e(b)
matrix V = e(V)

/* Storing them as variables */
capt drop xaxis beta se meanheterovar 
gen xaxis = .
gen beta = .
gen se = .
gen meanheterovar = .

forvalues i = 1/4{
	qui replace xaxis = `i' in `i'
	qui replace beta  = B[1,`i'] in `i'
	qui replace se    = sqrt(V[`i',`i']) in `i'
	qui sum age if heterovar_age_qt == `i' & year == $lastyearpre ///
	& bergen_commuters == 1
	qui replace meanheterovar = r(mean) in `i'
}

/* 95 % CI */
capt drop plusse minusse
qui gen plusse = beta + 1.96 * se
qui gen minusse = beta - 1.96 * se


/* === FIGURE ============================================================*/

twoway ///	
	(rcap plusse minusse xaxis ///
		, sort color(black)) ///			
	(scatter beta xaxis ///
		, sort lcolor(black) mcolor(black) lpattern(solid)) ///	
	, graphregion(color(white))	///
	xtitle("") ytitle("Treatment effect") ///
	name(fig3b, replace) ///
	xlabel(1 "Single w/o kids" 2 "Single w. kids" 3 "Couple w/o kids" 4 "Couple w. kids", angle(55)) ///
	legend(off) yline(0) ///
	ylabel(-0.02(0.02)0.08) /// 
	scale(1.3)

graph export "${figures}fig5b.png", as(png) replace
graph export "${figures}fig5b.pdf", as(pdf) replace



/* === Fourth heterogeneous treatment: education ================================== */		  

gen heterovar_educ = .
replace heterovar_educ = max_educ + 1 if year == 2014
bys famid_num: egen aux = mean(heterovar_educ)
replace heterovar_educ = aux
drop aux

capt drop bergen_commuters_post_*	
forvalues i = 1/5 {
	gen bergen_commuters_post_`i' = bergen_commuters_post * (heterovar_educ == `i')
}			  

eststo clear
reghdfe $yvar ///
	bergen_commuters_post_* /// DiDiD: post commuter in bergen
	i.year##i.heterovar_educ /// 
	$xvar ///
	i.commuters##i.heterovar_educ##i.year /// 
	i.bergen_commuters##i.heterovar_educ /// 
	i.bergen##i.heterovar_educ  ///
	, absorb( ///
		i.grk#i.year ///
		i.famid_num  ///
	) vce(cluster grk)	
eststo reg_3
forvalues i = 1(1)5 {
	summarize heterovar_educ if e(sample) == 1  & heterovar_educ ==`i'
	estadd scalar MeanHetero_`i' = r(mean), replace
	estadd scalar NHetero_`i' = r(N), replace
	summarize $yvar if e(sample) == 1  & year == $lastyearpre ///
	& bergen_commuters == 1 & heterovar_educ ==`i'
	estadd scalar MeanValue1_`i' = r(mean), replace
	summarize $yvar if e(sample) == 1  & year == $lastyearpost ///
	& bergen_commuters == 1 & heterovar_educ ==`i'
	estadd scalar MeanValue2_`i' = r(mean), replace
}
eststo reg_3 
estimates save ${ster}reg_hetero_${yvar}_educ, replace

/* === For making figures without running regression ======================== */
global reg2yvar cars
estimates use ${ster}reg_hetero_${yvar}_educ
/* ========================================================================== */

matrix B = e(b)
matrix V = e(V)

/* Storing them as variables */
capt drop xaxis beta se meanheterovar 
gen xaxis = .
gen beta = .
gen se = .
gen meanheterovar = .
forvalues i = 1/5{
	qui replace xaxis = `i' in `i'
	qui replace beta  = B[1,`i'] in `i'
	qui replace se    = sqrt(V[`i',`i']) in `i'
	qui sum heterovar_educ if heterovar_educ == `i' & year == $lastyearpre ///
	& bergen_commuters == 1
	qui replace meanheterovar = r(mean) in `i'
}

/* 95 % CI */
capt drop plusse minusse
qui gen plusse = beta + 1.96 * se
qui gen minusse = beta - 1.96 * se


/* === FIGURE ============================================================*/

twoway ///	
	(rcap plusse minusse xaxis ///
		, sort color(black)) ///			
	(scatter beta xaxis ///
		, sort lcolor(black) mcolor(black) lpattern(solid)) ///	
	, graphregion(color(white))	///
	xtitle("") ytitle("Treatment effect") ///
	name(fig3c, replace) ///
	xlabel(1 "Unknown" 2 "< high school" 3 "High school" 4 "Uni., short" 5 "Uni., long", angle(55)) ///
	legend(off) yline(0) ///
	ylabel(-0.02(0.02)0.08) /// 
	scale(1.3)

graph export "${figures}fig5c.png", as(png) replace
graph export "${figures}fig5c.pdf", as(pdf) replace



/* === Fifth heterogeneous treatment: work distance ================================== */		  

  
xtile heterovar_wd_qt = dist if year == 2014, nquantiles(5)
table heterovar_wd_qt, stat(max dist)
bys famid_num: egen aux = mean(heterovar_wd_qt)
replace heterovar_wd_qt = aux
drop aux

drop bergen_commuters_post_*	
forvalues i = 1/5 {
	gen bergen_commuters_post_`i' = bergen_commuters_post * (heterovar_wd_qt == `i')
}			  

eststo clear
reghdfe $yvar ///
	bergen_commuters_post_* /// DiDiD: post commuter in bergen
	i.year##i.heterovar_wd_qt ///
	$xvar ///
	i.commuters##i.heterovar_wd_qt##i.year /// 
	i.bergen_commuters##i.heterovar_wd_qt /// 
	i.bergen##i.heterovar_wd_qt ///
	, absorb( ///
		i.grk#i.year ///
		i.famid_num ///
	) vce(cluster grk)	
eststo reg_3
forvalues i = 1(1)5 {
	summarize dist if e(sample) == 1  & heterovar_wd_qt ==`i'
	estadd scalar MeanHetero_`i' = r(mean), replace
	estadd scalar NHetero_`i' = r(N), replace
	summarize $yvar if e(sample) == 1  & year == $lastyearpre ///
	& bergen_commuters == 1 & heterovar_wd_qt ==`i'
	estadd scalar MeanValue1_`i' = r(mean), replace
	summarize $yvar if e(sample) == 1  & year == $lastyearpost ///
	& bergen_commuters == 1 & heterovar_wd_qt ==`i'
	estadd scalar MeanValue2_`i' = r(mean), replace
}
eststo reg_3 
estimates save ${ster}reg_hetero_${yvar}_wd, replace
/* === FIGURE of the last of the regressions === */
/* Saving estimates */
matrix B = e(b)
matrix V = e(V)

/* Storing them as variables */
capt drop xaxis beta se meanheterovar 
gen xaxis = .
gen beta = .
gen se = .
gen meanheterovar = .
*local coef = $fig4lastyear - $fig4firstyear /* The number of coef we want to keep */
forvalues i = 1/5 {
	qui replace xaxis = `i' in `i'
	qui replace beta  = B[1,`i'] in `i'
	qui replace se    = sqrt(V[`i',`i']) in `i'
	qui sum dist if heterovar_wd_qt == `i' & year == $lastyearpre ///
	& bergen_commuters == 1
	qui replace meanheterovar = r(mean) in `i'
}

/* 95 % CI */
drop plusse minusse
qui gen plusse = beta + 1.96 * se
qui gen minusse = beta - 1.96 * se


/* === FIGURE ============================================================*/

twoway ///	
	(rcap plusse minusse xaxis ///
		, sort color(black)) ///			
	(scatter beta xaxis ///
		, sort lcolor(black) mcolor(black) lpattern(solid)) ///	
	, graphregion(color(white))	///
	xtitle("Work distance quintile") ytitle("Treatment effect") ///
	name(fig3e, replace) ///
	xlabel(1 "5-7 km" 2 "7-10 km" 3 "10-13 km" 4 "13-19 km" 5 "19-50 km", angle(55)) ///
	legend(off) yline(0) ///
	ylabel(-0.02(0.02)0.08) ///
	scale(1.3)

graph export "${figures}fig5e.png", as(png) replace
graph export "${figures}fig5e.pdf", as(pdf) replace

/* === Sixth heterogeneous treatment: public_diff_car_time ================================== */		  
 
xtile heterovar_ptdiff_qt = PublicDiffCarTime_fam_mean if year == 2014, nquantiles(5)
table heterovar_ptdiff_qt, stat(mean PublicDiffCarTime_fam_mean)
table heterovar_ptdiff_qt, stat(max PublicDiffCarTime_fam_mean)
bys famid_num: egen aux = mean(heterovar_ptdiff_qt)
replace heterovar_ptdiff_qt = aux
drop aux

drop bergen_commuters_post_*	
forvalues i = 1/5 {
	gen bergen_commuters_post_`i' = bergen_commuters_post * (heterovar_ptdiff_qt == `i')
}			  

eststo clear
reghdfe $yvar ///
	bergen_commuters_post_* /// DiDiD: post commuter in bergen
	i.year##i.heterovar_ptdiff_qt /// 
	$xvar ///
	i.commuters##i.heterovar_ptdiff_qt##i.year /// 
	i.bergen_commuters##i.heterovar_ptdiff_qt /// 
	i.bergen##i.heterovar_ptdiff_qt ///
	, absorb( ///
		i.grk#i.year ///
		i.famid_num ///
	) vce(cluster grk)	
eststo reg_3
forvalues i = 1(1)5 {
	summarize PublicDiffCarTime_fam_mean if e(sample) == 1  & heterovar_ptdiff_qt ==`i'
	estadd scalar MeanHetero_`i' = r(mean), replace
	estadd scalar NHetero_`i' = r(N), replace
	summarize $yvar if e(sample) == 1  & year == $lastyearpre ///
	& bergen_commuters == 1 & heterovar_ptdiff_qt ==`i'
	estadd scalar MeanValue1_`i' = r(mean), replace
	summarize $yvar if e(sample) == 1  & year == $lastyearpost ///
	& bergen_commuters == 1 & heterovar_ptdiff_qt ==`i'
	estadd scalar MeanValue2_`i' = r(mean), replace
}
eststo reg_3 
estimates save ${ster}reg_hetero_${yvar}_ptdiff, replace

/* === FIGURE of the last of the regressions === */
/* Saving estimates */
matrix B = e(b)
matrix V = e(V)

/* Storing them as variables */
drop xaxis beta se meanheterovar 
gen xaxis = .
gen beta = .
gen se = .
gen meanheterovar = .
forvalues i = 1/5 {
	qui replace xaxis = `i' in `i'
	qui replace beta  = B[1,`i'] in `i'
	qui replace se    = sqrt(V[`i',`i']) in `i'
	qui sum PublicDiffCarTime_fam_mean if heterovar_ptdiff_qt == `i' & year == $lastyearpre ///
	& bergen_commuters == 1
	qui replace meanheterovar = r(mean) in `i'
}

/* 95 % CI */
drop plusse minusse
qui gen plusse = beta + 1.96 * se
qui gen minusse = beta - 1.96 * se


/* === FIGURE ============================================================*/

twoway ///	
	(rcap plusse minusse xaxis ///
		, sort color(black)) ///			
	(scatter beta xaxis ///
		, sort lcolor(black) mcolor(black) lpattern(solid)) ///	
	, graphregion(color(white))	///
	xtitle("Quintiles of public transit quality") ytitle("Treatment effect") ///
	name(fig3f, replace) ///
	xlabel(1 "< 29 min" 2 "29-37 min" 3 "37-48 min" 4 "48-71 min" 5 "> 71 min", angle(55)) ///
	legend(off) yline(0) ///
	ylabel(-0.02(0.02)0.08) /// 
	scale(1.3)

graph export "${figures}fig5f.png", as(png) replace
graph export "${figures}fig5f.pdf", as(pdf) replace
